Binary black hole merger gravitational waves and recoil in the large mass ratio limit 
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Spectacular breakthroughs in numerical relativity now make it possible to compute spacetime 
dynamics in almost complete generality, allowing us to model the coalescence and merger of binary 
black holes with essentially no approximations. The primary limitation of these calculations is 
now computational. In particular, it is difficult to model systems with large mass ratio and large 
spins, since one must accurately resolve the multiple lengthscales which play a role in such systems. 
Perturbation theory can play an important role in extending the reach of computational modeling 
for binary systems. In this paper, we present first results of a code which allows us to model the 
gravitational waves generated by the inspiral, merger, and ringdown of a binary system in which one 
member of the binary is much more massive than the other. This allows us to accurately calibrate 
binary dynamics in the large mass ratio regime. We focus in this analysis on the recoil imparted to 
the merged remnant by these waves. We closely examine the "antikick," an anti-phase cancellation 
of the recoil arising from the plunge and ringdown waves, described in detail by Schnittman et al. 
We find that, for orbits aligned with the black hole spin, the antikick grows as a function of spin. 
The total recoil is smallest for prograde coalescence into a rapidly rotating black hole, and largest 
for retrograde coalescence. Amusingly, this completely reverses the predicted trend for kick versus 
spin from analyses that only include inspiral information. 
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I. INTRODUCTION AND BACKGROUND 



A. Modeling binary systems in general relativity 

After roughly three decades of effort, numerical rela- 
tivity can now model nearly arbitrary binary black hole 
configurations. Following Pretorius' pioneering "break- 
through" calculation [l], and then the successes of the 
Brownsville and Goddard groups using techniques that 
required only modest modifications to the methods they 
used before the breakthrough j2|, Q , the past few years 
have seen an explosion of activity. Recent work has stud- 
ied the impact of the many physicalparameters which 
describe binaries, such as mass ratio [J, [5|, spin and spin 
alignment and eccentricity [Tol.llll|. As numerical 

models have improved, analytic tools for modeling bi- 
nary systems [l2| and connecting numerics and analytics 
have likewise matured. In particular, the effective one- 
body (EOB) [l3''-^16'| approach to binary dynamics, which 
maps the dynamics of a binary to that of a point parti- 
cle moving in an "effective" spacetime corresponding to 
a deformed black hole, has been found to outstandingly 
describe the outcome of numerical relativity calculations 
after some adjustable parameters in the EOB framework 
are calibrated to numerical calculations [T7l - [20j . Our un- 
derstanding of the two-body problem in general relativity 
has never been better. 

These efforts are largely motivated by the need for 
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accurate models of coalescing black holes to detect and 
measure merger signals in the data of gravitational- wave 
(GW) detectors. Black holes with masses of roughly 
10^ — 10^ Mq indisputably reside at the cores of essen- 
tially every galaxy with a central bulge [2l|, . In the 
hierarchical growth of structure, these black holes will 
form binaries as their host galaxies merge and grow [23| ; 
estimates of how often such binaries form indicate that 
the proposed space-based detector LISA [23, HI] should 
be able to measure at least several and perhaps several 
hundred coalescences over a multiyear mission lifetime 
[26j . There is already a catalog of candidate binaries 
in this mass range, such as active galaxies with double 
cor es l27l - [29| , systems with doubly-peaked emission lines 
[sol ISlI). and systems that appear to be periodic or semi- 
periodic, such as the blazar OJ287 [S^j. The last year 
or so of the binary's life will generate GWs at frequen- 
cies to which LISA is sensitive; measuring those waves 
will make it possible to precisely map the distribution of 
cosmic black hole masses and spins, opening a new obser- 
vational window onto the high-redshift growth of cosmic 
structure. 

Less massive black hole binaries (several to several hun- 
dred Mq) will be targets for the ground-based GW de- 
tector network, currently including LIGO 33], Virgo [3^, 
and GEO [11], and hopefuUyincluding the proposed de- 
tectors LCGT [1^, AIGO and the "Einstein Tele- 
scope" [13 in the future. Formation scenarios and event 
rate estimates in this band are much less certain, since 
the demographics of the relevant black holes and scenar- 
ios for them to form binaries are not as well understood as 
in the supermassive range. However, scenarios involving 
dynamic binary formation in dense clusters suggest that 
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this network can plausibly expect an interesting event 
rate [39l - l43j , strongly motivating the construction of bi- 
nary merger models for these detectors. 

Finally, moving back to the LISA band, binaries in 
which one member is much less massive than the other 
are expected to be an important source. Such extreme 
mass ratio binaries are created when a stellar mass sec- 
ondary 1 — 100 Mq) is scattered through multibody 
interactions onto a highly relativistic orbit of a roughly 
10^ Mq black hole in the center of a galaxy. Though 
rare on a galaxy-by-galaxy basis, enough galaxies will 
be in the range of LISA that the number measurable is 
expected to be several dozen to several hundreds [il ]. 
The waves from these binaries largely probe the quies- 
cent spacetimes of their larger (presumably Kerr) black 
hole, making possible precision tests of the strong-field 
nature of black hole spacetimes [i^ . 

In short, astrophysical binary black holes will come 
in a wide range of mass ratios. Computational models 
must be able to handle systems with mass ratios ranging 
from near unity, to millions to one. Each mass m sets 
a lengthscale Gm/<? which the code must be able to re- 
solve. Large mass ratios require codes that can handle a 
large dynamic range of physically important lengthscales. 

Perturbation theory is an excellent tool for modeling 
binaries with very large mass ratios. In this limit, the 
binary's spacetime is nearly that of its largest member, 
with the smaller member acting to distort the metric 
from the (presumably) exact Kerr solution of that "back- 
ground." It is expected that tools based on perturbation 
theory will be crucial for modeling extreme mass ratio 
systems described above (mass ratios of 10^ : 1 or larger). 
Even for less extreme systems, perturbative approaches 
are likely to contribute important wisdom, working in 
concert with tools such as numerical relativity and the 
effective one-body approach. 

The foundational examples of such an analysis are 
the papers of Nagar, Damour, and Tartaglia [46] and of 
Damour and Nagar [13] ■ In that work, the EOB frame- 
work is used to construct the quasi-circular late inspiral 
and plunge of a small body into a non-rotating black hole. 
Regge- Wheeler- Zerilli methods [i^ 11^ are then used to 
compute the GWs that arise from a small body that fol- 
lows that trajectory into the larger black hole. Those 
authors use this large mass ratio system as a "clean lab- 
oratory" for investigating binary dynamics, and advocate 
using these techniques as a tool for probing delicate is- 
sues such as the form of the waves which arise from the 
plunge, and the matching of the final plunge waves to the 
late ringdown dynamics of the system's final black hole. 

Our goal here is to develop a similar toolkit based on 
perturbation theory applied to spinning black holes. We 
have developed two perturbation theory codes which we 
use to model different aspects of binary coalescence. Both 
codes solve the Teukolsky equation [50|, computing per- 
turbations to the curvature of a Kerr black hole. One 
code works in the frequency domain [5ll.[5^. which works 
well for computing the averaged flux of quantities such 



as energy and angular momentum carried by GWs. The 
other code works in the time domain (ssl . [s^ ] , which is ex- 
cellent for calculating the aperiodic GW signature of an 
evolving source. As originally proposed in Ref. (ssj . we 
have developed a hybrid approach which uses the best 
features of both the time- and frequency-domain codes 
to model the full coalescence process. (Although our ul- 
timate goal is to develop a set of tools similar to those 
developed by Damour, Nagar, and Tartaglia, we note 
that our techniques are the moment largely numerical, as 
opposed to the mixture of numerical and analytic tech- 
niques developed in Refs. 0,|43|. It would be worthwhile 
to connect the work we present here to the body of EOB 
work, but have not yet begun doing so in earnest.) 

As we were completing this paper, a perturbation- 
theory-based analysis of binary merger was presented by 
Lousto et al. [S^ . Their analysis does not use the Teukol- 
sky equation, but is otherwise very similar in style and 
results to what we do here. In particular, they note as 
we do here that the perturbation equations terminate the 
merger waveform in a set of ringdown waves in a very nat- 
ural way, thanks to the manner in which the equation's 
source redshifts away as the infalling body approaches 
the large black hole's event horizon. This behavior was 



also pointed out and exploited by Mino and Brink [57] 
in their (largely analytic) perturbative analysis of recoil 
from waves from the late plunge. We expand on this 
point in more detail at appropriate points later in the 
paper. 

As our use of the Teukolsky equation requires, we as- 
sume that a binary can be well-described by a small body 
moving in the spacetime of a (much larger) Kerr black 
hole. We first build the worldline that the smaller body 
follows as it slowly inspirals and then plunges into the 
black hole. We assume that, early in the coalescence, 
the small body moves on a geodesic of the background 
Kerr spacetime. Using the frequency-domain perturba- 
tion theory code to compute their rates of change, we 
allow the energy i?, angular momentum L^, and Carter 
constant Q of this configuration to evolve. (In fact, we 
confine ourselves to equatorial orbits in this analysis, so 
Q — Q throughout the binary's evolution.) This drives 
the smaller body in an adiabatic inspiral through a se- 
quence of orbits, until we approach the last stable orbit 
of the large black hole^. 

We then make a transition to a plunging orbit, us- 
ing the prescription of Sundararajan [s^j which in turn 
generalized earlier work by Ori and Thorne^ j59i] . By 



^ At present, we do not include the conservative impact of self 
forces. These forces are included in the EOB-framework analyses 
of Damour, Nagar, and Tartaglia. 

^ A similar approach to the transition from inspiral plunge, but 
valid for arbitrary mass ratios and presented using the EOB 
framework, was developed by Buonanno and Damour Il4( , and 
appeared in press before Ref. [59l |: we thank T. Damour for clari- 
fying this chronology to us. This approach is used in Refs. [4^ . |47|| 
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properly connecting the adiabatic inspiral to a plunge, 
we make a full worldline describing the small body's coa- 
lescence with the larger black hole. This worldline gives 
us the source for our time-domain perturbation theory 
code, from which we compute the GWs generated by the 
system as the small body evolves from the (initially near 
geodesic) inspiral through the plunge and merger. The 
waves which we compute in this way have qualitatively 
the same "inspiral, merger, ringdown" structure seen in 
numerical relativity simulations, though much work re- 
mains to quantify the degree of overlap. 

As an illustration of the utility of our perturbative 
toolkit, we focus in this paper on the problem of GW re- 
coil. Studies of GW recoil have been particularly active 
in recent years; we review this problem and its literature 
in the next subsection. 



B. Gravitational- wave recoil 

The asymmetric emission of GWs from a source carries 
linear momentum. The system then recoils to enforce 
global conservation of momentum. Early work demon- 
strated the principle of this phenomenon jsTI , l62j ; Beken- 
stein [m appears to have been the first to appreciate the 
important role it could play in astrophysical problems. 
Much recent work has focused on the recoil imparted to 
the merged remnant of binary black hole coalescence. 

The first estimates of binary black hole kick were made 
by Fitchett [6^] . He treated the gravitational interaction 
as Newtonian and included the lowest order mass and 
current multipoles needed for GW emission to compute 
the recoil velocity. This early calculation predicted that 
recoil velocities could approach thousands of km/s, which 
is greater than the escape velocity for many galaxies. Be- 
cause of his restriction to low-order radiation formulas, 
and his use of Newtonian gravity to describe binary dy- 
namics, it was clearly imperative that Fitchett's calcula- 
tions be revisited; a prescient analysis by Redmount and 
Rees [ill particularly argued for the need to account for 
the effect of black hole spins in the coalescence. 

Over the past several years, quite a few calculations 
have substantially improved our ability to model the re- 
coil in general relativity. The various approaches can be 
grouped as follows: 

• Black hole perturbation theory. As discussed exten- 
sively above, black hole perturbation theory is a 
good tool for describing binaries involving a mas- 
sive central black hole (of mass M) and a much 
less massive companion (of mass /x) . Shortly after 
Fitchett's pioneering binary calculation, Fitchett 
and Detweiler examined whether strong-field grav- 



to compute the transition from the slow, adiabatic inspiral to 
plunge. 



ity changed the conclusions using perturbation the- 
ory |66| . Twenty years later, Favata, Hughes, and 
Holz [63 argued that, properly extrapolated, rea- 
sonable results can be obtained for quantities such 
as the integrated black hole kick up to a mass ra- 
tio /i/Af ^ 0(0.1). Unfortunately, the Favata et 
al. analysis has a rather large final error since the 
frequency-domain tools they use do not work well 
at modeling the GWs arising from the final plunge 
of the smaller body into the large black hole. One 
of our goals in this analysis is to revisit that calcu- 
lation and reduce those substantial error bars. 

Another application of perturbation theory is the 
"close-limit approximation," [68] which describes 
the last stages of a merging binary as the dynamics 
of a distorted single black hole. Sopuerta, Yunes, 
and Laguna ^6£] applied the close-limit approxima- 
tion to describe the final waves from unequal mass 
binaries, obtaining results that compare very well 
with those that have since been computed within 
"full" numerical relativity. 

Finally, Mino and Brink [Ft'] used perturbative 
techniques to model the waves from the plunge, 
quantifying the manner in which the geometry of 
the final infall impacts the kick imparted to the 
binary. As already mentioned, their analysis also 
took advantage of the manner in which the source 
redshifts away as the infalling body approaches the 
larger black hole's event horizon. 

• Post-Newtonian (PN) theory. FN theory describes 
the spacetime and the motion of bodies in the 
spacetime as an expansion in the Newtonian grav- 
itational potential Gm/rc^ (where m is a charac- 
teristic system mass, and r a characteristic black 
hole separation). Blanchet, Qusailah, and Will [T^I 
used an approach based on this expansion to sub- 
stantially improve estimates from the recoil from 
the final plunge and merger; though consistent with 
the results from (63], they were able to reduce the 
error bars by a substantial factor. More recently, 
Le Tiec, Blanchet, and Will [7l| combined a post- 
Newtonian inspiral with a close-limit computation 
of the merger and ringdown to compute the re- 
coil for the coalescence of non-spinning black holes. 
This analysis is quite similar in spirit to the one we 
present here, though it does not use perturbation 
theory throughout. 

• Numerical relativity. Not long after it first became 
possible to model the coalescence of two black holes 
in numerical relativity, this became the technique 
of choice for computing black hole recoil. No other 
technique is well-suited to computing wave emis- 
sion and spacetime dynamics for very asymmetric, 
strong-field configurations which are likely to pro- 
duce strong GW recoils. Numerical relativity was 
needed to discover the so-called "superkick" con- 
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figuration: an aUgnment of spin and orbital angu- 
lar momentum which results in a kick of several 
thousand kilometers per second [72l - [73 |. In most 
configurations, the kick tends to be substantially 
smaller, peaking at a few hundred kilometers per 
second ^E^- 

• Effective one-body. As already described, EOB de- 
scribes a binary as a test body orbiting in the space- 
time of a "deformed" black hole, with the deforma- 
tion controlled by factors such as the mass ratio of 
the binary. Damour and Gopakumar fTS^ first ex- 
amined the issue how to compute recoil within the 
EOB framework, analytically identifying the major 
contributions to the recoil that accumulates over 
a coalescence, including the importance of the fi- 
nal merger and recoil waves in providing an "an- 
tikick" contribution. By calibrating some parame- 
ters of the EOB framework with results from nu- 
merical relativity, EOB has had great success gen- 
erating waveforms and recoil velocities that match 
well with those from numerical relativity [t^, H^I • 

With the exception of the "superkick" configuration, all 
of these techniques predict recoils that peak at roughly a 
few hundred kilometers per second (depending on mass 
ratio, spins, and spin orbit orientation; see |8l| for de- 
tailed discussion and statistical analysis). This is sub- 
stantially lower than the peak predicted by Fitchett's 
original calculation; his overestimate can be ascribed to 
neglect of important curved spacetime radiation emission 
and propagation effects. 

In addition to their potential astrophysical applica- 
tions, recoil computations serve another important pur- 
pose: They are a common point of comparison for these 
four approaches to strong field gravity. The recoil veloc- 
ity from a merging binary is calculated by integrating the 
emitted radiation over some number of orbits. Any sig- 
nificant systematic error in the approach used used will 
tend to magnify the error in the estimated recoil veloc- 
ity. Thus, the evaluated recoils for a range of BH spins 
and mass ratios serve as a good platform for comparing 
various approaches to strong field binary models. 

C. This paper 

Our goal is to revisit and improve the estimate of black 
hole recoil via black hole perturbation theory that was 
originally developed in Ref. [s^]- That analysis predicts 
upper and lower bounds which are rather widely sepa- 
rated. This is because the analysis of could not ac- 
curately model wave emission from the final plunge and 
merger. Using the time-domain perturbation theory code 
developed and presented in Refs. [H, [13], we can now 
compute the contribution of those waves. As we describe 
in more detail in Sec.|Vl doing so completely reverses the 
conclusions of Ref. [67] regarding how the kick behaves 
as a function of spin. In particular, including the plunge 



and merger is crucial to correctly computing the "an- 
tikick," the out-of-phase contribution to the recoil that 
arises from the merger's final GWs. This contribution to 
a binary's total recoil was first identified and character- 
ized by Schnittman et al. [82,]. We find that the inability 
to include this contribution in Ref. [g^] is largely respon- 
sible for the large error bars in that analysis. 

We begin by reviewing in Sec. |ll] how we construct 
the worldline which the smaller member of our binary 
follows as it spirals into the larger black hole. As briefiy 
described above, we break this trajectory into a slowly 
evolving "inspiral" (Sec. Ill A[) followed by a transitional 
regime (Sec. Ill B|) that takes the binary into a final plunge 
and merger (Sec. Ill C|) . This review is left general, so that 
in principle one could describe this dynamics for generic 
orbital geometry. We specialize in our analysis here to 
the simplest circular and equatorial orbits (Sec. HID]) . 

We next briefiy review how we compute gravitational 
radiation from a body moving on this trajectory. As 
mentioned above, our approach is based on finding solu- 
tions to the Teukolsky equation fssl for Kerr black hole 
perturbations. We review this equation's general prop- 
erties in Sec. IIIII and then discuss the principles behind 
solving it in the frequency domain (Sec. lIII"A|) and in the 
time domain (Sec. IIITB|) . Section llVI summarizes how 
one computes the radiation's linear momentum and the 
recoil of a merged system. 

Section |V] presents the results of our analysis. We be- 
gin in Sec. IV Al with general considerations on how our 
results scale with mass ratio. Because we work strictly 
within the context of linearized perturbation theory, all 
of our results can be easily scaled to different mass ra- 
tios, provided that the scaling does not change the sys- 
tem so much that the validity of perturbation theory 
breaks down. Reference argued that a modified scal- 
ing would allow us to estimate with reasonable accuracy 
quantities related to the recoil even out of the perturba- 
tive regime. Although those arguments are valid during 
the adiabatic inspiral, they break down when the mem- 
bers of the binary merge. 

In Sec. IV Bl we then discuss in some detail the grav- 
itational waveform we find for binary coalescence in the 
large mass ratio limit. We examine the different mul- 
tipolar contributions to the last several dozen cycles of 
inspiral, followed by the plunge and merger. These ex- 
amples illustrate the manner in which the coalescence 
waves very naturally evolve into a "ringdown" form. As 
discussed in some detail in Sec. lIIIBl this behavior arises 
by virtue of how the Teukolsky equation's source term 
goes to zero, so that its solutions transition to their ho- 
mogeneous form, as the infalling body approaches the 
large black hole's event horizon. Mino and Brink [s^] 
first appear to have exploited this behavior, which was 
also seen in recent work by Lousto, Nakano, Zlochower, 
and Campanelli This demonstrates the power of 

perturbative methods at modeling physically important 
aspects of the merger waves. 

Section IV CI examines the recoil that arises from these 
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waves, focusing on how it depends (for the circular, equa- 
torial case that we study) on the spin of the larger black 
hole. This analysis demonstrates very clearly the im pac t 
of the "antikick" first reported by Schnittman et al. [S^] ■ 
With the antikick taken into account, the smallest re- 
coils come from the largest spins when the merger is a 
prograde sense; the largest spins come from retrograde 
mergers with large spins. The waves which give the sys- 
tem its antikick come from those produced by the final 
plunge and merger, demonstrating very clearly the sub- 
stantial impact these waves have on the system. We con- 
clude this section by briefly discussing the convergence 
of our recoil results as a function of black hole spin. In- 
terestingly, we find that the number of modes we must 
include in order for our results to converge is a strong 
function of the black hole's spin — rapid spin, prograde 
cases need more modes than do slow spin cases, which in 
turn need more modes than rapid spin, retrograde cases. 
We conclude the paper by discussing in Sec. lVII how these 
tools may be used to expand the reach of two-body mod- 
eling in general relativity, and our future plans. 

Throughout our analysis we generally use units in 
which G — c— 1. We sometimes use c = 3 x 10^ km/sec 
in order to present kicks in "physical" units. 



II. BUILDING THE INSPIRAL AND PLUNGE 
TRAJECTORY 



Roughly speaking, our coalescence model has two in- 
gredients. First, we compute the worldline that the 
small body follows as it spirals from large radius through 
plunge into the black hole. We then use that worldline to 
build the source for the Teukolsky equation and compute 
the GWs that are generated as the smaller body follows 
the worldline into the black hole. Though for simplic- 
ity we describe these ingredients as though they stand in 
isolation, they are in fact strongly coupled. We describe 
here how we compute the inspiral and plunge trajectory, 
deferring discussion of how we compute radiation from 
this trajectory to Sec. IIIII Throughout, we indicate how 
these steps are coupled to one another. 

The trajectory which the small body follows can be 
broken into three pieces: An early time inspiral, in which 
the smaller member of the binary is approximated as 
evolving through a sequence of bound orbits of the larger 
black hole; a late time plunge, in which the small body 
falls into the larger black hole; and an intermediate tran- 
sition which smoothly connects these two regimes. We 
now briefly review how we model these different pieces. 

In all of these regimes, we treat the zeroth order mo- 
tion of the small body as a geodesic of the Kerr spacetime. 
These geodesies must be augmented by the conservative 
action of a self force if one's goal is to make a model that 
faithfully reproduces the phase of binary black hole GWs. 
For our present goal of estimating the GW recoil, we ex- 
pect that the error due to neglecting this force is not im- 
portant. Kerr black hole geodesies [85J are described by 



the following equations for the motion in Boyer-Lindquist 
coordinates r, 6, (j), and t: 

(2.1) 
(2.2) 
(2.3) 
(2.4) 

The potentials appearing here are 

R = [E{a^ + r'^) -aL^]^ 

-A [(L, - aS)2 + + Q] , (2.5) 
Ve = Q- cos^ 6 [a^ (^^ _ _^ ^,^^2 g,^2j ^ ^^^.Q) 

= csc^ BLy^ -aE^^{E (r^ + a^) - L^a] (2.7) 
Vt ^ a {L^ - aE sin^ 0) 
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The quantity M is the large black hole's mass, a is that 
hole's Kerr spin parameter, and fj, is the mass of the 
smaller body which perturbs the black hole spacetime. 
The functions S = -|- cos^ 6 and A = - 2Mr + 
a^. In the absence of radiation emission, the energy E, 
axial angular momentum L^, and Carter constant Q are 
constants of the motion; up to initial conditions, choosing 
these three constants defines a geodesic. 

Equations (|2.1[) - (|2.4p are the starting point for build- 
ing the smaller body's inspiral and plunge worldline. We 
now describe in some detail how we use them for this 
computation. 



A. The inspiral 

We approximate the inspiral as a slowly evolving se- 
quence of bound Kerr geodesies (neglecting for now con- 
servative aspects of the self interaction). Momentarily 
ignore the impact of radiation emission. In this limit, the 
orbits are determined by selecting E, L^, and Q plus ini- 
tial conditions, and are completely characterized by three 
orbital frequencies describing their periodic motions in 
the r, 9, and (j) coordinates [11]. This periodic nature 
means that functions built from the orbital motion can 
be usefully represented by a discrete Fourier expansion. 

To build our inspiral, we assume that radiation acts 
slowly enough that, to a good approximation, we can 
treat the small body's worldline as a Kerr geodesic at 
each moment. We then use the frequency-domain Teukol- 
sky solver described in Sec. IIIII to compute the rates at 
which E, Lz, and Q evolve due to GW backreaction. 
From these rates of change, we build the time-varying 
parameters E{t), Lz{t), and Q{t) which describes the se- 
quence of orbits the small body passes through on its 
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inspiral. More detailed discussion of this procedure is 
given in Refs. (E^ and @. 



B. The last stable orbit and the transition to 
plunge 

Our assumptions, and hence our procedure for com- 
puting the inspiral, break down as the small body ap- 
proaches the last stable orbit, or LSO. This is worth de- 
scribing in some detail. For bound orbits, the function 
R{r) defined in Eq. (|2.5|) generally has four real roots. 
Denote these roots ri > r2 > > r4. The root is 
generally inside the event horizon'^, and is not interest- 
ing for our discussion. The roots ri, r2, and r^, on the 
other hand, are quite important. When these roots are 
distinct, the geodesic describes an eccentric orbit that os- 
cillates between ri (apoapsis) and r2 (periapsis). When 
ri = r2 > 7*3, the geodesic describes a circular orbit at 
r = ri. (In this case, we also have dR/dr = at r = ri.) 
When r2 = r^, the orbit is marginally stable. (The triple 
root ri = r2 = denotes a marginally stable circular 
orbit.) Once we reach this point, the small body will 
rapidly plunge into the black hole. This condition de- 
fines the LSO. 

As inspiral proceeds, the roots r2 and r^ approach 
one another, indicating that GW backreaction is car- 
rying the small body toward the LSO. We model the 
transition from slowly evolving geodesies through the 
LSO to plunge by expanding the equations of motion 
around their behavior at the LSO, as described in Ref. 
[ssj (which generalizes Ref. [l^). More specifically, we 
take the constants in the transition to be given by 



E{t) ~ -ElSO + (i ^ ^LSo)£'LSO 



(2.9) 



and similarly for Lz and Q. Here, -Elso and -Elso are 
the energy and its rate of change at the LSO (the latter 
calculated using our frequency-domain Teukolsky equa- 
tion solver), and ^lso is the time at which the LSO is 
reached. We integrate the geodesic equations using this 
form from a time tstart < ^lso until a time tend > ^lso- 
Reference [11] describes how we choose istart and t^nd as 
a function of parameters such as the black hole spin a 
and binary mass ratio. For our purposes, it is enough to 
note that, provided they are chosen within a well-defined 
range, our results are robust to that choice — varying 
istart and tend docs not significantly change the recoil. A 
more careful investigation may clarify the optimal way 
to define these transition parameters. 



C. The plunge 

For t > ti^so, the geodesies described by E{t), Lz{t), 
and Q{t) correspond to plunging geodesies, i.e., trajecto- 
ries which fall into the large black hole. As described in 
Ref. [Ill, the transition matches onto a plunging trajec- 
tory most simply by just holding these parameters con- 
stant for t > tend- This is justified by the fact that ra- 
diation reaction does not have a strong impact in the 
final plunge [1^ 113, : careful analysis indicates that 
an orbit's energy and angular momentum remain nearly 
constant during the final plunge into the black hole. 

As the small body approaches the black hole, its mo- 
tion as viewed by distant observers appears to "freeze" 
onto the generators of the event horizon**. When this 
happens, the source term of the Teukolsky equation red- 
shifts to zero [cf. Eq. (2.46) of Ref. [53]]. Since the ho- 
mogeneous Teukolsky equation's solutions are the quasi- 
normal modes of the binary's large black hole, this means 
that the final cycles of radiation from our coalescing sys- 
tem are very naturally given by the system's ringdown 
modes. 



D. Specialization to circular equatorial orbits 

Until now, we have kept the discussion of these inspi- 
ral and plunge trajectories general in order to emphasize 
that our approach can be applied to totally generic co- 
alescences. For this first analysis, we now focus on the 
simplest interesting case, circular orbits confined to the 
equatorial plane of the larger black hole. In this limit, 
geodesic orbi_ts_are totally characterized by the orbit's 
radius; 



gives an outstanding summary of their 



Ref. 
properties. 

The energy, angular momentum, and Carter constant 
of circular equatorial orbits are given by 



E = 



1 - 2M/r ± aM^/^/r^^^ 
v/1 -3Af/r±2aMi/2/r3/2' ' 



(2.10) 



Q 



± / ' ' — ^ , 2.11 

^1 -3M/r±2aMV2/r3/2 

. (2.12) 



Upper sign refers to prograde coalescences (orbital an- 
gular momentum parallel to black hole spin), lower sign 
to retrograde (antiparallel). These orbits are character- 
ized by a single frequency associated with the azimuthal 
motion. 



•3/2 ± aMi/2 • 



(2.13) 



^ In fact, r4 = for equatorial orbits {Q = 0) and for orbits of 
Schwarzschild black holes (a = 0). 



** This behavior led many researchers to call these solutions "frozen 
stars" in the early literature. 
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FIG. 1: Example quasi-circular inspiral and plunge trajectory. For this calculation, the larger black hole's spin was set to 
a — 0.3A/, and the binary's mass ratio is ^/M = 10"''. Left panel shows r{t), the trajectory's radius as a function of time; 
right panel shows the same trajectory as viewed in the equatorial plane of the larger black hole. On the left, the dotted line 
at r = 4.98Af labels the radius of the prograde last stable orbit. Our trajectory, which starts at r = 5.23M, executes about 
25 orbits before crossing this point. The inset there zooms in on the region t ~ 2260M, showing the smaller body's trajectory 
as it approaches the event horizon at r = 1.955M. On the right, the heavy circle at r = 1.955A/ is the hole's event horizon; 
notice how the particle quickly "locks" onto the horizon after the plunge which follows its slow inspiral. 



The last stable orbit is located at 

nso/M 



3 + ^2 T V(3-^i)(3 + ^1+2^2) , 

(2.14) 

(l + a/M)i/3 + (i_a/M)i/3l ^ (2.15) 



Zo = 



(2.16) 



Our procedure for building an inspiral and plunge tra- 
jectory reduces, in the equatorial circular limit, to the 
following algorithm: 

1. Choose a mass /i for the smaller body, and mass 
M and spin a for the larger black hole^. Pick an 
initial orbital radius r and an initial azimuth (j) for 
the smaller body's trajectory. 

2. Evolve through a sequence of circular, equato- 
rial orbits using E computed with the frequency- 
domain code (described in the following section). 
For these orbits, we don't need to compute Q since 



^ As is common with codes of this sort, we normalize most dimen- 
sionful quantities to M. As such, we really pick mass ratio /^/A/; 
the impact of M can then be accounted for in post-processing 
after the numerics have been evaluated. 



Q — Q over the entire sequence. Also, in this case 
Lz is simply related to E, so computing it does not 
provide additional information. 

3. As we approach the last stable orbit, switch to the 
transition trajectory following Ref. |58{ . In partic- 
ular, Ref. [58] describes how to choose the times at 
which we start and end the transition regime, which 
depends in detail on the system's mass ratio, the 
larger black hole's spin, and the orbit geometry. For 
circular, equatorial orbits, this prescription reduces 
to that given in Ref. [s^. 

4. When we reach icnd, hold the parameters E and 
Lz constant, and allow the small body to follow 
the plunging trajectory so defined into the larger 
black hole. 

An example trajectory is shown in Fig. [1] For this fig- 
ure, we examine a binary with a mass ratio n/M — 10"'*. 
The larger black hole has a spin a = 0.3Af. We start our 
prograde trajectory at r = 5.23M, close to the LSO at 
r = 4.98Af. The smaller body orbits roughly 25 times 
before crossing the LSO; shortly thereafter, it rapidly 
plunges into the black hole, locking onto the horizon as 
seen by distant observers. The inset in the left panel of 
Fig. [1] zooms in on its approach to the horizon, showing 
that our plunge trajectory smoothly asymptotes to the 
final "horizon locking" behavior. 
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III. COMPUTING RADIATION 

We compute gravitational radiation from our model 
binaries using the Teukolsky equation, whicli describes 
the evolution of curvature perturbations to a Kerr black 
hole [il,!!!]. In Boyer-Lindquist coordinates, it is given 
by 



: svcL I 



dtt'^ —dt4,'^ 



-2s 



la cos 
1 



smf 



A 



1 



a 
A 



sm 

- (s^ cot^e- s) * 



a{r - M) 
A 



-47r fr^ + cos^ I 



sin2 6' 



(3.1) 



where M is the mass of the black hole, a its angu- 
lar momentum per unit mass, A = — 2Mr + o?, 
r± = M±VAf2 — a^. The quantity s is the "spin weight" 
of the field under study. Choosing s — Q means the per- 
turbing field is a scalar field; s — ±.1 describes a spin-1 
(electromagnetic) perturbation, and s — ±2 describes 
gravitational perturbations. For s = +2, the field \E' 
is given by the Weyl curvature scalar ■00 (see [s^ for 
precise definitions and discussion of this quantity); for 
s = —2, — {r — ia cos 9)'^tfj4, where V'4 is another cur- 
vature scalar. We use s = —2 since ip4 is a natural choice 
to study outgoing radiation. Once ■04 is known, we then 
know the GWs that the binary produces, since 



1p4 



1 /d^h. 



(3.2) 



far from the black hole. 

The T on the right hand side of Eq. p.ip is a source 
term constructed from the stress-energy tensor describing 
a point-like body moving in the Kerr spacetime. This 
stress-energy tensor is given by 



lap 



laUf} 



(3.3) 



= ^^^f^6[r-r{t)]6[e-e{t)]6[^ 
Y.t sm 



(3.4) 

On the top line, x denotes an arbitrary spacetime event, 
z(t) describes the worldline that the point-body follows, 
and r is proper time along that worldline. On the second 
line, we have performed the integral and written the re- 
sult in terms of the body's motion in the Boyer-Lindquist 
coordinates r, 9, and 0, parameterized by coordinate time 
t. On both lines, Ua = dza/dr is the 4- velocity of the 
body as it moves along its worldline. 

Note in particular the t = dt/dr that appears in the de- 
nominator of Eq. p.4p . This is the time-like component 
of the smaller body's geodesic motion, described by Eq. 



p.4p . As the small body approaches the horizon, t — > oo 
— the passage of coordinate time (time as measured by 
distant observers) diverges per unit proper time as mea- 
sured by that body. This is the mechanism by which the 
source term "redshifts away" as the small body falls into 
the large black hole, smoothly converting the Teukolsky 
equation into its homogeneous form. 

The source T is constructed from this Tap by project- 
ing onto a tetrad that describes radiation, and then ap- 
plying a particular integro-differential operator; see Refs. 
[52, 53, 83] for detailed discussion of its nature. For our 
purposes here, the key thing to note is that we must 
construct the worldline which the small body follows in 
order to compute the radiation associated with its mo- 
tion. Different approximations are appropriate to dif- 
ferent regimes of the coalescence, which is why we have 
developed two rather different codes for solving Eq. (13.11) . 
We now briefiy summarize the techniques behind these 
two codes, and how we use our solutions. 



A. Radiation in the frequency domain 

As was originally found by Teukolsky [83], Eq. (|3.1 
separates. For s — —2, we put 



{ra(f)—ujt) 



^4 = 7 — ^ [ dujy^ Ri,nujir)Shni9)e'' 

(r — la cos Or J ^-^ 

hn 

(3.5) 

The function Sijn{9) is a spin- weighted spheroidal har- 
monic, and can be constructed by expanding on a basis 
of spin- weighted spherical harmonics f5l| . The function 
Rimu {r) is found by solving a second-order ordinary dif- 
ferential equation. Its limiting behavior is 

i?im^(r^oo)ocZ°°e^"'^* , (3.6) 

corresponding to purely outgoing radiation far away, and 

i?/™^(r^r+) cxZ^e-^'^'-* , (3.7) 

corresponding to purely ingoing radiation on the event 
horizon. The wavenumber k = co — mu}+, where a;+ = 
a/2Mr+ is the angular velocity of the hole's event hori- 
zon. In both of these equations, r* is the so-called "tor- 
toise coordinate," 



2Mr^ 



■In 



2M 



2Afr_ 



■In 



2M 



(3.8) 

The ingoing and outgoing solution is thus characterized 
by the coefficients Z°° and . For details of how we 
compute these numbers, see Refs. [sills^. 

This frequency-domain approach to solving Eq. p.ip 
is most useful when the source T has a discrete frequency 
spectrum. The function ■04 can then be written as a sum 
over harmonics of the source's fundamental frequencies. 
This is the case for geodesic black hole orbits; see Ref. 
(52i] for extensive discussion. 
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For the circular, equatorial case, orbits and hence the 
source T are completely characterized by the frequency 
defined in Eq. (|2?T3l) . The frequency uj in Eq. (|331) 
becomes mJl^, and the coefficients Z°° and are then 
determined by uj and the harmonic indices / and m. Once 
those coefficients are known, it is not difficult to compute 
the rates of change of E, L^, and Q. The coefficients 
Z°°'^ are labeled by the indices I and m, and we have 



E 



hn 



(3.9) 
(3.10) 



where w„j — mQ^. Strictly speaking, the I sum appear- 
ing here is from Z = 2 to infinity, and m is from —I to I; 
in practice, the sums converge to double precision accu- 
racy once I is of order a few to a few dozen, depending 
on how fast the smaller body orbits. See Refs. [5l|, [52| 
for extensive discussion of convergence issues, as well as 
for discussion of how to compute the down-horizon con- 
tribution to the rates of change. Also, see Ref. [9^ for 
discussion of how to compute the rate of change of Q. 



the time domain with less than a 1% error compared to 
a frequency-domain code over a large span of orbital pa- 
rameter space. 

We use the transition and plunge trajectory described 
in the previous section to provide the worldline z(t) 
and 4- velocity Ua- As we have already highlighted, the 
Teukolsky equation (|3.ip becomes homogeneous at late 
times thanks to the manner in which t — oo as the in- 
falling body approaches the event horizon. When T = 0, 
the solutions of Eq. (|3.ip are the larger black hole's quasi- 
normal modes modes. This means that the late-time so- 
lution in our coalescence model is dominated by quasi- 
normal modes of the larger black hole. By virtue of aris- 
ing in a natural way from the behavior of our source 
term, these modes are properly phase connected to the 
preceding inspiral and plunge waves. 



IV. COMPUTING RECOIL FROM RADIATION 

Once we have computed ip4, it is not difficult to com- 
pute the rate at which linear momentum is carried by 
the waves. Letting T^^p^ denote the Isaacson [9lj stress- 
energy tensor for GWs, we have 



B. Radiation in the time domain 

Because they work best when the source has a discrete 
frequency spectrum, we only use frequency-domain tech- 
niques to describe the inspiral, when the system is accu- 
rately described as slowly evolving through a sequence of 
orbits. When this description is not accurate (such as in 
the final plunge, or when inspiral is sufficiently rapid that 
the system does not spend many cycles near a given or- 
bit), Eq. (|3.5p is ill-suited to describing solutions of the 
Teukolsky equation. To handle this case, we solve Eq. 
p.ip directly in the time domain. 

In the code we have developed for this, we take advan- 
tage of the Kerr spacetime's axial symmetry to write the 
field ^' as [H m 

^{t,r,e,cj)) ^ ^e™V^„,(i,r,0) . (3.11) 

m 

Equation p.ip is then solved as a (2+l)-dimensional par- 
tial differential equation for the modes $m- 

The major difficulty in numerically solving Eq. p.ip is 
coming up with a good description of the source term. 
One challenge is to represent a point-like source on a nu- 
merical grid [we use finite difference techniques to solve 
Eq. p.ip ]. In Refs. [s^ and [U, we have developed a 
discrete representation of a delta function which works 
very well on a finite-difference grid. This function is de- 
fined so that our representation of the delta function and 
of its first two derivatives preserves various integral iden- 
tities. For cases in which a comparison can be made 
(e.g., for non-evolving generic geodesic orbits), we find 
that this representation allows us to compute GWs in 



dP\t) 
dt 



lim / n' T^"^ dVL 

r— >-oo 



= lim — 



= lim 



dh+ 
~dt 

fdh+ 
\ dt 
fdh^ 



dhx 
dt 



dfl 



.dhy. 



lim — 

r^oo 4:7T 



\ dt 



ipidt' 



.dhx 



dn 



dn 



(4.1) 



The quantity denotes the Cartesian direction vector 
in the large radius limit. 



n" 



sm W cos ( 
sin 9 sin (, 
cos 9 . 



(4.2) 
(4.3) 
(4.4) 



This quantity must then be integrated over time to find 
the momentum carried by the GWs: 



P\t) 



dP^t')^. 



dt 



dt' . 



(4.5) 



Imposing global conservation of momentum, the recoil 
velocity of the system is then given by 



<,(i) = -P\t)/M 



(4.6) 



Equation ()4.6|) will be our primary tool for computing 
black hole kicks in this analysis. 
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For the inspiral, which we model using frequency- 
domain methods, these formulas reduce to fairly simple 
results. In the r — >■ cx) limit, 

^4 = J E ^^nSim (0)e™'(^-"**) . (4.7) 

Inserting this expansion, the momentum flux formula 
(|4.1[) reduces to a sum over overlap integrals between dif- 
ferent modes of the radiation field. This integral sharply 
constrains the mode numbers which contribute to this 
formulas sum. For Schwarzschild, the integral over 9 can 
be expressed as a Clebsch-Gordan coefficient, and we find 
Z' G [I — 1,1^1 -\- 1]; a similar but more complicated result 
describes the integral for Kerr. For any spin, we find 
to' = m ± 1 for P^(t) and P^(t) [P^(i) = for the equa- 
torial orbits we consider here] . Details of this calculation 
will be presented in a separate analysis [92|. For the fi- 
nal plunge and merger portions of the coalescence, we 
simply evaluate Eqs. (|4.ip . (j4.5p . and (14.61) using the ^^4 
computed with our time-domain code. 

V. RESULTS: THE COALESCENCE 
WAVEFORM AND RECOIL 

We now put the pieces of this formalism together to 
compute the waveforms from binary black hole coales- 
cence and to calculate recoil. We begin by describing 
some issues with extrapolating from the truly perturba- 
tive mass ratios we study here (Sec. IV Al) . and then de- 
scribe in more detail how we assemble the full inspiral 
trajectory and its associated waveform (Sec. IV Bp before 
discussing our results for the recoil fSec. IIV[) . As already 
mentioned, we focus in this analysis on quasi-circular 
equatorial configurations. We conclude (Sec. IV Pp by 
describing the convergence of our recoil results. We find 
that as we go to large spin, prograde mergers may require 
a large number of m modes [cf. the axial decomposition 
p. lip we use] to give convergent results. 

A. Mass ratio dependence considerations 

By using the Teukolsky equation to model coalescence, 
we are by construction working to first order in mass 
ratio — the curvature scalar tp^ that we compute ne- 
glects all corrections of order (/x/M)^. Since the various 
fluxes we compute (energy, momentum, angular momen- 
tum) follow from the modulus squared of "04, it likewise 
follows that these fluxes are all strictly proportional to 
{ji/MY . The recoil velocity, as an integral of the momen- 
tum flux, should likewise scale essentially with (/x/M)^. 
We may expect small deviations from this scaling since 
the timescales of inspiral and of plunge and merger do 
not scale with mass ratio in quite the same way. How- 
ever, we have found that a (fj./M)'^ scaling describes the 
final recoil very accurately for all mass ratios more ex- 
treme than /i/M = 10~^; we have not examined mass 



ratios less extreme than this yet. Since the scaling with 
mass ratio is trivial, we will present detailed results for 
only one choice, ^i/M = 10^^. In our summary figure 
for total recoil velocity as a function of spin (Fig. [S]), we 
normalize our results by the scaling (/i/M)^. 

In Ref. [H^l, it was argued that one can improve the 
ability of perturbation theory to extrapolate out of the 
perturbative regime by replacing the (/i/il/)^ which de- 
scribes the momentum flux and the recoil velocity with 

In this argument, it is claimed that in extrapolating out 
of the perturbative regime it is useful to interpret the 
small body's mass /i as the system's reduced mass, and 
the large black hole's mass M as the system's total mass. 
A similar interpretation of these masses has been shown 
to give excellent results interpreting the head-on colli- 
sions of black holes [88j. The function /(/i/M) has a 
maximum /max = 0.01789 at ^/M = 0.2 (corresponding, 
after remapping the meaning of these mass parameters, 

to TOsmall/Wiarge = 0.382). 

As we will discuss in more detail later in this section, 
using this scaling doesn't work quite as well as we might 
have hoped. The key issue is that in our perturbative 
framework, we assume there exists a stationary back- 
ground spacetime which we can expand around, and that 
this background does not evolve during the coalescence. 
This means, for example, that a binary which contains a 
large Schwarzschild black hole at early times will evolve 
to a single Schwarzschild black hole at late times; we fail 
to account for the evolution of this black hole's spin dur- 
ing the merger. This is a minor error when the mass ratio 
is small, but is significant for large mass ratio. In partic- 
ular, for mass ratios fJ./M 0.1 or larger, the spin of the 
final black hole will change substantially in the merger. 
By not evolving the spin properly, we do not get the late 
time spectrum of merger /ringdown waves correct, with 
important consequences for the system's final kick. 

B. Example waveform 

Figures [2l|3l and |4] present coalescence waveforms for a 
binary with /i/M = 10"**, and in which the larger black 
hole has spin a/M = 0.6, 0, and —0.6 respectively. We 
focus on the late waves, including the final plunge and 
ringdown. The data for these figures were generated us- 
ing the time-domain code discussed in Sec. IIIIBI In the 
largest panel, we show the wave including all contribu- 
tions with |to| < 6; the three smaller panels show indi- 
vidual contributions from the to = 1, 2, and 3 modes. 

In all three examples we show, the general character 
of the waveforms is essentially the same: a slowly evolv- 
ing chirping sinusoid that terminates in an exponentially 
damped ringdown. Two aspects of the waveforms clearly 
differ as we move from a — 0.6M to a = — 0.6M. First, 
notice that the waveform for a — 0.6M is clearly of rather 
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FIG. 2: Coalescence waveform computed with perturbation 
theory for a binary with mass ratio fJ./M = 10~*, and in 
which the larger black hole has spin a = O.GAf. We show 
the + polarization of the waveform as viewed in the binary's 
equatorial plane; the x waveform is zero from this viewing 
angle. The wave is normalized by D, the distance from source 
to observer, and the origin of the time axis is arbitrary. The 
waveform shown in the top panel includes contributions from 
all modes with \m\ < 6; the individual contributions for m = 
1, m = 2, and m = 3 are shown below. Notice how the 
smoothly chirping "inspiral" waves blend naturally into the 
rapidly damped "ringdown" which terminates this waveform. 



higher frequency than for a = 0, which in turn is higher 
than for a = —0.6M. This is not surprising, and is a 
simple consequence of the orbit's geometry: the LSO, 
which approximately dehneates the transition from in- 
spiral to final plunge, is at tlso = 3.83Af for a = 0.6M, 
rLso = 6M for a = 0, and tlso = 7.85Af for a = -0.6M. 
As the LSO moves to larger radius, the orbital frequency 
associated with it sweeps lower. 

Second, the final ringdown waves damp more quickly 
as we move from the prograde to the retrograde configu- 
ration. This is also not surprising, and follows naturally 
from the damping behavior of a Kerr black hole's quasi- 
normal modes: modes which are "parallel" to a hole's 
spin (i.e., have m > for a > 0, and vice versa) are 
much more long lived than "antiparallel" modes. See, 
for example. Fig. 45 of Ref. [89] (noting that Chandra's 
sign convention on the Fourier transform means that 

"^Chandra = " ™us)- 



FIG. 3: Same as Fig. (2] but the large black hole has spin 
a = in this case. The frequencies which describe this wave 
are generically lower than those shown in Fig. [2] since the 
transition from inspiral to plunge happens at larger radius 
thanks to smaller spin parameter in this binary. In addition, 
the final ringdown waves damp out more rapidly than in the 
a — 0.6M case. 
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FIG. 4: Same as Figs. [2] and [S] but now for black hole spin 
a = — 0.6M (i.e., same black hole as in Fig. [21 but now for 
a retrograde orbit geometry). The frequencies characterizing 
this wave are again lower than for the two previous examples, 
and the ringdown waves damp even more rapidly. 



C. Recoil versus spin 

Figure[S]summarizes how the kick imparted to a binary 
behaves as a function of spin for mass ratio n/M = 10"**. 
In this plot, we show the magnitude of the recoil that has 
accumulated up to some time t. Since the origin of the 
time axis is not particularly interesting, we have shifted 



the various tracks so that we can easily compare how the 
recoil varies as a function of spin. 

The clearest feature apparent here is that, especially 
for prograde coalescences (a > 0), the recoil grows to 
some large positive value, but then is strongly suppressed 
by an "antikick" to something significantly smaller. The 
suppression is a very strong function of the large black 
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a = 0.6M, 1.5 for a = 0.3M, 1.2 for a = 0, and 1.02 for 
a — — 0.3M. The late-time kick shown in Fig. [5]is nicely 
fit by the formula 

w,ec(a)/c ~ [0.0440 - 0.0099(a/M) - 0.0114(a/Af )2 
-0.0312(a/M)3] ifi/Mf . (5.7) 

Over most of the relevant parameter space, this comes in 
right between the "upper" and "lower" estimates of Ref. 
[H [compare to Eqs. (1) and (2) of Ref. [93]]- 

The antikick behavior we see agrees at least qualita- 
tively with the trends seen in Schnittman et al. 182']. It 
is hard to calibrate the quantitative agreement between 
these analyses, since (as discussed above in Sec. IV Ap 
extrapolating from the perturbative regime into that of 
the mass ratios considered in Ref. is not as sim- 
ple as the arguments in Ref. [67| would suggest. Con- 
sider the Schwarzschild coalescence results. If we use the 
[ti/Mf -> fifJ./M) rule suggested in Ref. we find 



FIG. 5: Summary of recoil versus time for large black hole 
spins a/M € [-0.9, -0.6, -0.3, 0, 0.3, 0.6, 0.9]. For each track, 
the time origin is arbitrary, so we have shifted the data in 
order to cleanly display all seven recoil trends shown here. 
The key feature we find is the manner in which (especially for 
large spin, prograde coalescences) the kick builds to a large 
positive value, followed by an "antikick" that brings the to- 
tal accumulated recoil down to much smaller values. The 
antikick is especially strong when the spin is large and the co- 
alescence is prograde, and is essentially non-existent for large 
spin retrograde coalescence. 



hole's spin: For the five cases which show antikick behav- 
ior in Fig. m the peak kick tiP^'^'^ and the late time kick 
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0.027 in/M)^; (5.3) 



/c = 0.039 ifi/M)^ ; (5.4) 



(5.5) 



(a = 0) = 0.044/(^/M)c 
< 235 km/sec . 



(5.8) 



On the second line, we have used /max = 0.01789 in order 
to estimate how large the kick can be in this case. 

Taken at face value, this suggests that the recoil of 
Schwarzschild black holes has a maximum of 235 km/sec, 
34% higher than the maximum value of 175 km/sec found 
in careful numerical relativity calculations ^94] . However, 
in those numerical relativity calculations, the final black 
hole is not Schwarzschild, but has a spin a ~ 0.67M. 
Figure [S] tells us we should expect a larger "antikick" 
when the final black hole is rapidly spinning. Applying 
the same extrapolation to our recoil data for a = 0.6M 
(the nearest value to a = 0.67M in our dataset) leads to 



..late 



(a = 0.6M) = 0.027/(^/A/)c 
< 144 km/sec . 



(5.9) 



This is about 18% lower than the numerical relativity 
prediction. The lesson we take from this is that naive 
extrapolation from the small mass ratio regime does not 
give a good estimate of the final kick in the comparable 
mass case. Because the background spacetime is fixed, 
we do not accurately describe the system's final state and 
hence the last waves that it emits during the coalescence. 



-0.3M: 



yPeak 
.,,latc 



/c 



0.048 in/Mf 
0.047 ifi/M)^ 



(5.6) 



In other words, we find that the antikick suppresses the 
maximum recoil by a factor of 53 for a = 0.9Af , 2.7 for 



D. Convergence 



As discussed in Sec. IIIIB| our time-domain perturba- 
tion theory code expands the field \E' in axial modes; cf. 
Eq. (|3.1ip . The angular integral in Eq. (|4.ip takes the 
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form 



dt 



dt 



2tt 



J2 I d<lycosq^e'^">^^^,^t„ 



) (5.10) 



oc 



J2 I d0sin0e'(™-"')^$„ci>;^, 



X! ('^(™+l),m' - ^(m-l),m') ^m*m' (^-H) 



Contributions from terms with m = m' vanish; the recoil 
arises from beating between adjacent m- modes. 

The question we now address is how many m-modes 
must be included in order to accurately compute the re- 
coil. We have found that this is a strong function of black 
hole spin: When the black hole has large positive spin, 
many more modes are needed for convergence than for 
small or retrograde coalescence. 

Tables HI - IIIII summarize convergence data for three of 
the cases presented in Fig. [5] We show, as a function 
"T-max [the value of m and m' at which the sums in Eqs. 
(|5.10p and (|5.1ip are terminated] the peak magnitude of 
the momentum flux normalized by (/i/Af)^, 



V 



y/{dP^/dt)^ + {dPy/dty 



(5.12) 



The "max" subscript means that we select the maximum 
of this quantity over the timespan for which we compute 
the momentum flux. This quantity is given in units of 
M~^. We also show the percentage change in V as we 
increase Wmax by one. 



TABLE I: Convergence of recoil with azimuthal mode for 
a = — 0.6M. First column is mmax, the largest value of m 
we include. Second column is the value of V, the peak mag- 
nitude of the momentum flux normalized by (fi/M)^. The 
third column gives the percentage change in V we find as we 
increase m,nax from the previous value. 





r 


% change 


2 


2.855 X 10" 


-3 




3 


4.030 X 10" 


-3 


29.2% 


4 


4.557 X 10" 


'3 


11.6% 


5 


4.807 X 10" 


-3 


5.2% 


6 


4.930 X 10" 


-3 


2.5% 



Tables |T] - IIIII indicate that, once several modes have 
been computed, the fractional error in the momentum 
flux decreases by roughly a factor of two with each unit 
increase in m. However, the magnitude of the relative 
error is a rather strong function of black hole spin. For 
a — — 0.6M, we find that going from mi„ax = 5 to 
"T-max — 6 changes the momentum flux by only 2.5%. In- 
cluding additional modes presumably will only produce 



TABLE II: Convergence of momentum flux with m for a = 0. 
All details are as in Tab. IH 





P 


% change 


2 


1.712 X 10"^ 




3 


4.188 X 10"^ 


58.9% 


4 


5.508 X 10"^ 


24.0% 


5 


6.182 X 10"^ 


10.9% 


6 


6.532 X lO-^' 


5.4% 



TABLE III: Convergence of momentum flux with m for a 
0.6 A/. All details are as in Tab. HI 



^^max 


P 


% change 


2 


1.373 X 10" 


'3 




3 


7.488 X 10" 


-3 


81.7% 


4 


1.105 X 10" 


-2 


32.2% 


5 


1.302 X 10" 


^2 


15.1% 


6 


1.412 X 10" 


-2 


7.8% 



percent-level changes. For a = 0.6M by contrast, the 
flux changes by nearly 8% as mmax is increased from 5 to 
6. Many modes are clearly needed to accurately compute 
the waves (and the recoil from these waves) as the large 
black hole's spin approaches the Kerr maximum. 



VI. CONCLUSIONS AND FUTURE WORK 

Now that numerical relativity has effectively solved the 
two-body problem in general relativity, a major task for 
researchers has become to explore the parameter space 
of binary coalescence. This will insure that wave mod- 
els constructed as templates for GW data analysis fully 
encompass the range of behaviors that are likely in real 
binary mergers, and allow us to more fully understand 
the phenomenology of binary black hole merger astro- 
physics. In this analysis, we have demonstrated that per- 
turbation theoretical techniques based on the Teukolsky 
equation are an excellent tool for extending the reach of 
our computations, allowing us to model large mass ra- 
tios that are challenging for 3-1-1 numerical simulations, 
but may be of astrophysical significance. Our analysis 
joins previous work by Damour and colleagues |46|. ml. 
Mino and Brink 



521 



and by Lousto and colleagues [56| 
which likewise used perturbation theory to model large 
mass ratio binaries. By using the Teukolsky equation, we 
can explore how the larger black hole's spin impacts the 
analysis, exemplified by our demonstration of how the 
previously identified "antikick" [s^ strongly depends on 
this spin. 

Two directions for future analysis strike us as partic- 
ularly noteworthy. First, the major motivation for this 
work is that perturbation theory makes exploring param- 
eter space computationally fast and simple. As such, it 
would be worthwhile to continue this exploration, exam- 
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ining how the waveform varies as a function of spin-orbit 
alignment, and exploring (for example) how the antikick 
evolves as one varies the inclination smoothly from the 
prograde to the retrograde geometry. Preliminary calcu- 
lations of this behavior indicate that the antikick rapidly 
evolves with spin-orbit alignment, consistent with the re- 
sults of Mino and Brink 15 which demonstrate a strong 
dependence on the final kick with the plunge geometry. 

Second, as Damour and Nagar have emphasized [43], 
particularly useful application comes by including input 
from the effective one-body formalism in our description 
of the small body's motion; input from perturbation the- 
ory can likewise be used to calibrate certain parameters 
in the EOB framework. Now that the spin-au gme nted 
Hamiltonian for binary systems is understood [95l , [96l |. 
we expect that work to extend EOB to more broadly in- 
clude the impact of spin will become very active. We 
hope that the tools we have presented here will be use- 
ful for further refining what has already proved to be a 
valuable tool for modeling coalescing binaries. 
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